% [T1Hat, bHat, aHat, residual] = rdNlssasha(data, nlsS)
%
% Finds estimates of T1, a, and b using a reduced-dimension nonlinear least
% squares approach. The model a+b*exp(-t/T1) is used, where a and b are
% complex and T1 is real.
% The residual is the rms error between the data and the fit.
%
% INPUT:
% data - the data to estimate from
% nlsS - struct containing the NLS search parameters and
%        the data model to use
%
% written by J. Barral, M. Etezadi-Amoli, E. Gudmundson, and N. Stikov, 2009
%  (c) Board of Trustees, Leland Stanford Junior University

function [T1Est, bEst, aEst, res,sd,r2] =  SASHA(data, nlsS,fitnum)
try
    [T1Init bInit aInit res] = rdNlsPr(data,nlsS);
catch
    T1Init = 1000;
    aInit = max(data);
    bInit = -aInit;
end
[tVec,order] = sort(nlsS.tVec);
data = squeeze(data);
data = data(order);
data = data(:);
polluted_pos = zeros(1,length(data));
for ix = 2:length(data)
    if(data(ix)<= data(ix-1))
        polluted_pos(1,ix)  = 1;
        polluted_pos(1,ix-1)  = -1;
    end
end
%initially fitting
ti = tVec( 1:end );
dat = data(  1:end  );
aInit =  max(data);
T1Init = 1400;
if(fitnum == 2)
    x0 = [ aInit  T1Init];
else
    x0 = [ aInit  1 T1Init];
end
[aEst, bEst,T1Est]= subfitting(ti,dat,x0);
modelValue = aEst*(1- bEst*exp( - ti/T1Est ));
res = 1/sqrt(length(ti))*norm(1 - modelValue./dat(:));
R2=corrcoef(modelValue(:),dat(:));
if(abs(R2(2,1))<0.95)
    T1Est = 0;
end

% 
% if( sum(abs( polluted_pos ))> 0)
%     [Binary,Pos] = find( polluted_pos==-1 );
%     
%     for ix = 1:length( Pos)
%         ixPos = Pos(ix);
%         res = zeros(1,2);
%         for i = 1:2
%             ixPos = ixPos + ( i -1);
%             
%             temp_Polluted_pos = polluted_pos;
%             temp_Polluted_pos(1,ixPos) = 0;
%             
%             ti = tVec( temp_Polluted_pos ==0 );
%             dat = data( temp_Polluted_pos ==0 );
%             [aEst, bEst,T1Est]= subfitting(ti,dat,x0);
%             modelValue = aEst + bEst*exp( - ti/T1Est );
%             res(1,i) = 1/sqrt(length(ti))*norm(1 - modelValue./dat(:));
%         end
%         %make decision to remove which value
%         if(res(1,1)>= res(1,2)  )
%             polluted_pos(1,ixPos) = 0;
%         else
%             polluted_pos(1,ixPos-1) = 0;
%         end
%         %fit fit with the second data
%     end
% end
% 
% ti = tVec( polluted_pos ==0 );
% dat = data( polluted_pos ==0 );
% [aEst, bEst,T1Est]= subfitting(ti,dat,x0);
% 
% modelValue = aEst + bEst*exp( - tVec/T1Est );
% res = 1/sqrt(nlsS.N)*norm( (modelValue- data)./max(data));


%r2 = (sum(modelValue.*data(:)) -1./nlsS.N*sum(modelValue).*sum(data(:))).^2./( (sum(modelValue.^2) -1./nlsS.N*sum(modelValue).^2).*(sum(data(:).^2) -1./nlsS.N*sum(data(:)).^2));
sdi = median(abs( modelValue - data))/0.6745 ;
s=corrcoef([modelValue,data]);
r2=s(1,2);  %%������������е����ϵ��

% Delt_y_a = @(A,B,td,T1)sum( 1);
% Delt_y_b = @(A,B,td,T1)sum( exp(-td/T1));
% Delt_y_t1 = @(A,B,td,T1)sum(  B*exp(-td/T1)*td*(1/T1^2) );
% Di = zeros(3,3);
% for i = 1 :length(data)
%     dya = Delt_y_a(aEst,-bEst,tVec(i),T1Est );
%     dyb = Delt_y_b(aEst,-bEst,tVec(i),T1Est );
%     dyt1 = Delt_y_t1(aEst,-bEst,tVec(i),T1Est );
%     Di = Di + [dyt1*dyt1,dyt1*dya dyt1*dyb; ...
%         dyt1*dya,dya*dya dya*dyb; ...
%         dyt1*dyb,dya*dyb dyb*dyb] ;
% end
% D_half = Di./(sdi^2);
% C = inv(D_half );
sd =0 ; %sqrt(C(1,1));

    function [aEst, bEst,T1Est]=subfitting(ti, data,x0 )
        if( length(x0) == 3 )
            try
                x = fminsearch( ...
                    @(x)sum( abs( data- (   x(1)*(1-x(2)*exp(-ti/x(3))))).^2 ), ...
                    x0,optimset('display','off'));
                
            catch
                disp('there is an error');
                x=x0 ;
            end
            aEst = x(1);
            bEst = x(2);
            T1Est = x(3);
        else
            try
                
                x = fminsearch( ...
                    @(x)sum( abs( data- (   x(1)-x(1)*exp(-ti/x(2)))).^2 ), ...
                    x0,optimset('display','off'));
            catch
                disp('there is an error');
                x0
                x = x0;
                x
            end
            aEst = x(1);
            bEst = -x(1);
            T1Est = x(2);
        end
    end



% function f=model(x)
% f=x(1).*(1-exp(-ti/x(2)));
% end
%
% function f=lsfunction(para)
% f=sum((model(para)-dat).^2);
% end

end




